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The tracking of space objects requires frequent and accurate monitoring for collision 
avoidance. As even collision events with very low probability are important, accurate pre- 
diction of collisions require the representation of the full probability density function (PDF) 
of the random orbit state. Through representing the full PDF of the orbit state for orbit 
maintenance and collision avoidance, we can take advantage of the statistical information 
present in the heavy tailed distributions, more accurately representing the orbit states with 
low probability. The classical methods of orbit determination (i.e. Kalman Filter and its 
derivatives) provide state estimates based on only the second moments of the state and 
measurement errors that are captured by assuming a Gaussian distribution. Although the 
measurement errors can be accurately assumed to have a Gaussian distribution, errors with 
a non-Gaussian distribution could arise during propagation between observations. More- 
over, unmodeled dynamics in the orbit model could introduce non-Gaussian errors into 
the process noise. A Particle Filter (PF) is proposed as a nonlinear filtering technique 
that is capable of propagating and estimating a more complete representation of the state 
distribution as an accurate approximation of a full PDF. The PF uses Monte Carlo runs 
to generate particles that approximate the full PDF representation. The PF is applied in 
the estimation and propagation of a highly eccentric orbit and the results are compared 
to the Extended Kalman Filter and Splitting Gaussian Mixture algorithms to demonstrate 
its proficiency. 


I. Introduction 

In history, the pioneers of astronomy have always been curious about the planets by demonstrating 
attempts in characterizing their trajectories. Astronomers such as Johannes Kepler and Carl Friedrich Gauss, 
were among the first in the derivation of the laws of planetary motion and the methods of preliminary orbit 
determination respectively. 1 In order to accurately describe the trajectory of the orbit over time, we need 
some observations of the motion of the space object, good models of the orbit’s mathematical dynamics and 
measurements and an effective filtering method. 

One of the earliest filtering methods developed was a linear filtering method known as the Kalman Filter 
(KF), developed by Rudolph Kalman, as a linear filtering method for estimating the state by predicting its 
mean and covariance through the linear dynamics and updating them whenever measurements are available. 2 
The major assumption for this method was that the underlying mathematical model dynamics were linear 
and that all the error terms and measurements possess a Gaussian distribution. However, since most systems 
are nonlinear in nature, a derivative of the KF was derived, known as the Extended Kalman Filter (EKF) 3 to 
tackle nonlinear systems. This is probably the most widely used estimation algorithm for nonlinear systems 
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based on the assumption that the systems possess almost linear properties based on the time scale of the 
updates. 4 ~' The EKF works by computing the linearized dynamics and measurement models based on the 
assumption that the Jacobian matrix exists. a The mean of the state is predicted through the nonlinear 
dynamics and the covariance matrix is predicted linearly based on the computed state transition matrix. 

Other filtering methods that have been proposed for nonlinear dynamics include the Unscented Kalman 
Filter (UKF), 5,10 which uses a deterministic sampling technique called the unscented transform to select 
a minimal set of sample points known as the sigma points around the mean of the state. These sigma 
points are then propagated through the nonlinear dynamics and then the mean and the covariance are 
eventually constructed from these weighted sigma points. 5,11 The UKF is considered to have superior 
implementation properties to the EKF since it does not require the computation of the Jacobian matrices 
and it is comparable to the second order Gauss filter. 11, 12 However the UKF is still somewhat limited for non 
Gaussian distributions. Despite that the sigma points are propagated through the nonlinear dynamics, the 
final mean and covariance are calculated from these points and hence represent only second order moment 
statistics that are insufficient for a non Gaussian representation in some cases. Other variations of the KF are 
the square root unscented Kalman Filter, 12 central difference Kalman Filter, 13 ensemble Kalman Filter, 14 
Schmidt-Kalman Filter 15 and the Invariant Extended Kalman Filter 16 , to name a few. These filters attempt 
to estimate the state and the covariance by approximating the nonlinear dynamics with some representation 
of the mean and covariance while assuming Gaussian characterization. 

In order to incorporate the non-Gaussian or nonlinear elements in our models, we consider using the 
Particle Filter (PF). The PF is a sequential Bayesian algorithm that can be used for nonlinear systems and 
non Gaussian distributed states. The Bayesian approach to state estimation involves the construction of the 
PDF of the current state of an evolving system conditioned on the accumulated time history of observations or 
measurements. 17, 18 The central idea of the PF is to estimate the required state probability density function 
(PDF) by using a set of random samples (particles) with their corresponding weights. As the number of 
particles increases, the representation of the required PDF becomes more accurate. The propagation and 
prediction rules are such that the combined weights of the particles in a given region will approximate the 
integral of the PDF over the region to a value of one. 1 '’ 19 

Some work has been done by other researchers for orbit determination using the PF such as Gang and 
Xiao-Jun 20 by adapting the number of particles required, due to the fact that the PF experiences a problem 
known as sampling degeneracy. This is when the effective number of particles required to fully capture the 
PDF decreases over longer durations of time, caused by small measurement noise or data outliers. Their 
work was done to compare the PF to EKF and UKF over shorter durations of time and to deduce to optimal 
number of particles required for an orbit with a low eccentricity of 0.001165. In our study, we have found 
that orbits with lower eccentricities are well behaved over relatively shorter durations compared to highly 
eccentric orbits whose nonlinear dynamics are more pronounced at periapse. Hence, we attempt to address 
these effects in our work. A thorough investigation of the PF was also presented by Arulampalan et ah, 19 to 
study the various derivatives of the PF algorithm. Their analysis was focused solely on a nonlinear dynamic 
scalar case and over relatively shorter durations of time as well. In, 21 Soto also presented some work on a self 
adaptive PF that calculates the number of particles required at each iteration and illustrates an example for 
a 2-dinrension dynamic model. However, such constraints are challenging to apply to a larger dimensional 
model and are computationally costly. In relation to orbit determination, other work has also been done 
using the PF by Hwang and Speyer 22 on an adaptive resampling method of avoiding sampling degeneracy 
for relative positioning using GPS carrier-phase measurements. The orbit in question was nearly circular 
and hence do not depict huge nonlinear behavior effects that are present at apoapse and periapse. Also, the 
simulations were performed over relatively shorter durations. 

The research to be presented in this paper, is focused on applying the PF for statistical orbit determination 
by tackling highly nonlinear orbits with larger eccentricities and larger state dimensions. We demonstrate 
the state predictions and measurement update performances comparing the results obtained from using the 
PF to the EKF and a GMM that uses the UKF as a filter. 

II. Orbit Dynamics and Measurement Models 

Practical application of the orbit determination (OD) process requires the knowledge of the dynamic 
model and the observation model as accurately as possible. This includes deriving the equations of motion 

a since some systems may contain discontinuities or singularities, 8, 9 so the Jacobian may not always exist. 
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of the satellite, including the best known representations of many complex perturbations (e.g. atmospheric 
drag, solar pressure and higher-order terms in the gravity potential). However, since our purpose is to study 
and demonstrate the PF for this preliminary stage of research, we will work with an idealized assumption 
that only point mass gravitational forces exist between the two bodies (the earth and the satellite). 1,23 
Process noise will be included to account for un-modeled perturbations. The derivation for the governing 
equations are shown as follows: 


V(E0 

(1) 


(2) 

+ 1 j + Z k) 

(3) 


where U is the potential energy for the earth point mass and V is the mathematical operator for the gradient, 
defined as V = + jfd + ^k. The gradient of the potential energy equals to the force applied to the mass 

of the satellite. Hence, in deriving the equations of motion of the satellite, the acceleration of the satellite 
X is equated to the gradient of the potential energy, noting that p is the earth’s gravitational constant and 
the distance R of the satellite from the earth is defined as R = \fX 2 + Y 2 + Z 2 . The state vector X is the 
position and velocity states in Earth Centered Inertial (ECI) frame X = [IY Z V x V y V z ]. 

The observation model uses range measurements obtained from either an earth-based or space-based 
observing instrument within its field of view. The measurement used in this study is the range measurement, 
defined as the distance between an earth or space-based instrument performing the measurement and the 
satellite. If the position vector of the instrument is rj and the position vector of the satellite is r, the range 
Pj is defined as the scalar magnitude of the position vector of the satellite with respect to the instrument 
j, 23 given as 


Pj = [( r — r j) • ( r — r j)] ^ +w(t) (4) 

It is important that the correct transformations are implemented so that the instruments and the satellites 
observed are in the same frame. We present the dynamic model in simplified canonical form and the 
observation model as nonlinear expressions that are functions of the time varying state X. The process noise 
v(t) and measurement noise w(t) are added to the dynamic and observation models f and h respectively to 
model the uncertainties as zero-mean Gaussian random variables with covariances Q (t) and R(£). Equations 
(3) and (4) , can be written in the canonical form, where 

X = f (X, t) + v(f) (5) 

Y = h(X,t)+w(t) (6) 

Equations (5) and (6) are then presented as discrete models showing the evolution of the state sequence 
and measurements at discrete time intervals. 


Xfc = ffc(Xfc_i, Vfc_i) (7) 

Y fc = h fe (X fc ,w fc ) (8) 


III. Extended Kalman Filter 


The EKF is implemented based on the estimation of the states using the discrete system in Equations 
(7) and (8). These equations are linearized locally about the previous estimate as a sufficient description 
of the nonlinear system with the major assumption that the linearized dynamics are Gaussian distributed, 
where F and H are the Jacobian matrices for the nonlinear functions f and h respectively 


dfh (X) 

(9) 

dh k (X) 

IX — Xfc|fc_i 

(10) 


(11) 
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The EKF can be described in three steps; Initialization, Time Update or Prediction and Measurement 
Update, with detailed equations shown in. 4,19 The initialization step assumes the knowledge of the initial 
state condition with an a priori covariance matrix. The time update predicts the state and covariance 
matrix at the next time step and finally the state vector and covariance matrix are updated using the 
available measurements at that same time step. 

The EKF however, is not an optimal estimator model, this can cause the filter to diverge quickly. As only 
the mean and the state covariance matrix are propagated, the statistical description of the state estimate is 
limited to a Gaussian PDF. 

IV. Splitting Gaussian Mixtures Method 

The Splitting Gaussian Mixtures Method uses an initial Gaussian distribution to represent the state and 
splits this into a sum of Gaussian distributions upon detection of nonlinearity in the model during state 
propagations. This algorithm was developed in detail by DeMars, Bishop, Jah, Giza and Kelecy in 24-26 and 
was proposed to capture nonlinearities during state propagations and is also capable of using a separate filter 
such as the UKF for measurement updates. The PDF of the Gaussian Mixture Model (GMM) is defined as 
a sum of Gaussians as follows: 

L 

p(X) = mi.Pi) (12) 

i = 1 

p g (X: m,P) = (27r)"t|p|-| exp{-^(X - m) T P^ 1 (X - m)} (13) 

where m and P are the mean and covariance of the random vector X, L represents the number of components 
in the GMM and uy are the weights of the individual Gaussian components. The sum of all the weights need 
to add up to one to retain the support of the PDF property that the area underneath the function must 
equal one. 

The nonlinearity during the state prediction is detected by calculating the differential entropy that is 
equivalent to the amount of disorder or surprisal in a random variable or vector. Whenever the differential 
entropy exceeds a preset threshold, then nonlinearity has been detected during the state propagation. 24 The 
differential entropy is defined by 

H(X) = K{— logp(X)} (14) 

upon further evaluation and simplification, the differential entropy becomes 

tf(X) = ilog|27reP| (15) 

The time rate change of the differential entropy is given as 

H(X) = ^trace{p- 1 P} (16) 

where 

P(f) = F(X(i), t)P(t) + P(f)F T (X(f), t) (17) 

After further simplification, and on the special assumption that the linearized dynamical system has the 

property that the trace of F(X(t),t) = 0, then the differential entropy is a constant 

H{X) = 0 (18) 

and so it is easy to see that over time, when a nonlinearity is detected the assumption of the trace of 
F(X(f),f) will not hold and so the algorithm will trigger a splitting process to split the initial Gaussian 
distribution into a sum of Gaussians. The differential entropy for the linear case is computed using the 
covariance matrix obtained from computing the Jacobian and the state transition matrices. To compare 
this to the nonlinear case, the covariance matrix can be computed using the UKF. Moreover, there are cases 
when the F(X(f), t) ^ 0 and so both the linear and the nonlinear differential entropies must be computed at 
each time step online and this can be computationally costly. Moreover, the GMM method represents a part 
of the filtering algorithm. For this example the GMM uses the Unscented Kalman Filter (UKF) (please see 
References 5, 10 for details of the UKF algorithm) to perform the measurement update process which requires 
a number of sigma-points based on the dimension of the state for each individual Gaussian component. 
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V. Particle Filter 


In order to incorporate the higher order moments of the state and measurement errors, we need to be 
able to study the evolution of the full state PDF in the OD process. The Particle Filter (PF) is effective in 
estimating the full state PDF of nonlinear models or non-Gaussian PDFs. 19,20,2 ' The PF is a simulation- 
based filter based on the sequential Monte Carlo approach. The central idea of a PF is to represent the 

(i) N 

required probability density function (PDF) by a set of N » 1 random samples (particles) {X}. }. =1 , with 
associated weights such that 


N 

p(X k \Y k ) - X< <} ) (19) 

i- 1 

where the particles are identically and independently distributed (i.i.d) samples drawn from the initial 
PDF with an initial Gaussian distribution assumption. The weights are the probability values that are the 
approximations of the relative densities of the particles such that their sum equals one. 19 ' 20,2 ' The state 
estimate is given as a conditional density called the posterior density p(X k \Y k ). Estimators based on this 
posterior density are constructed from Bayes’ theorem, which is defined as 


P(X|Y) - (20, 

where p(X) is called the prior density (before measurement), p(Y|X)is called the likelihood and p(Y)is 
called the evidence (normalizes the posterior to assure the integral sums up to unity). 19,2 ' Let X k = 
(X 0 , X 1? ..., Xfc) and y k = (Yo,Y u ...Y k ) be the stacked vectors of states and observations up to time k 
and {w k ,i = 1,2,..., TV} represent the weights of the N particles at each time k. In the execution of the 
general PF algorithm, a common effect ensues where the variance of the weights increases over time. This 
phenomenon is known as sampling degeneracy, in that there is an insufficient number of effective particles 
required to fully capture the PDF. A way of overcoming the sampling degeneracy is to derive the importance 
distribution that minimizes the variance of the distribution of the weights. The corresponding minimum 
variance weights are given by: 


tu fc _ip(Y fe |X fe _i) 

(21) 


(22) 


This means that the minimum variance weights can be calculated before the particles are propagated 
to time k. From the equations above, we see that there is a problem in the sense that the evaluation of 
the integral generally does not have an analytic form and that the particles weight update are based on the 
knowledge of the state at the past time step and the current measurement only i.e. p(X k \X k -i,Y k ). So, a 
reasonable choice for the importance distribution is the transition prior p(Xfc|Xfc_i). 19 

The estimation process can be sectioned into three parts: initialization, time update or prediction and 
the measurement update. Each step is described below for the PF implementation. 

1. Initialization at k = 0, for i = 1,...,N 

• Sample N particles from the prior Xq ~ p(X g) 

• Calculate the weights of the particles from the initial distribution (assume Gaussian) w l 0 = 
P(Y„|X* 0 ) 

• Normalize the weights Wg = ^ using total weight wt = w l 0 

2. Time Update or Prediction at k > 1, for i = 1,...,N 

• Predict the particles through the dynamics 

Xfc = l) + v k (23) 

p(X fe |Y fc -i) = J p(X fe |X fe _ 1 )p(X fe _ 1 |y fc _i)dX fe _ 1 (24) 
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3. Measurement Update 

• Update the weights using the innovation from the measurements assuming that the measurements 
are Gaussian distributed 


w k = w k _ x - 


p( Y fc |Xt)p(X[.|Xj._ 1 ) 

q (x.i\xLM 


(25) 


• We assume that the evidence (also known as the importance distribution) q(X. l k \X k _\, y k ) in this 
case is equal to the prior density p(X. k \X. k _ 1 ) , 27 so that, 


w i k =w i k _ 1 xp(Y k \xi) 


(26) 


• Normalize the weights w l k = ^ using total weight wt = w k 

• The final posterior density estimate p(Xfc|Yfc) = Yf- T =1 w^S(X k -X^) 


Due to the fact that the prior density is usually broader than the likelihood, only a few particles will 
be assigned higher weight during the measurement update step causing sampling degeneracy. The solution 
to overcome sampling degeneracy is to resample the particles. Let jV e // denote the effective sample size, 


which can be estimated by: 19,28 N e ff = 




Also let N t h denote a lower threshold for the effective 


sample size, which is arbitrarily chosen with respect to the accuracy desired; the larger the threshold, the 
more accurate the PDF results and the more computational cost incurred. 

If N e ff > N t h, then the sampling degeneracy is not low enough, the filter continues on to step 2 for the 
next time update. Otherwise, we replicate the particle with the highest weight to replace the particles falling 
below the threshold and then normalize the weights. However, resampling does have its own shortcomings, 
since the particles that have higher weights are statistically selected each time thus losing the diversity of the 
samples. This loss of diversity may also cause the divergence of the estimates. To avoid the loss of diversity 
the replicated particles are “jittered” by adding process noise to spread the resampled particles . 19,28 


VI. Simulation Scenario 

For the first test case scenario in comparing the EKF to the PF, we consider a highly eccentric orbit of 
eccentricity 0.8181 with an orbital period of 24hrs. The satellite has four opportunities for measurements 
from 2 Deep Space Network ground stations located at Canberra, Australia and Madrid, Spain for a duration 
of 75 minutes each using the Satellite Toolkit (STK) software for a simulated realistic scenario. The last two 
measurement opportunities are obtained from a space-based observation network called TDRSS (Tracking 
and Data Relay Satellite System) for a duration of 15 minutes each. 


TDRS-X and TDRS-3 



Epoch: 04-Nov-2011 16:00:00 
Period: 24 hrs 


DSN: Canberra, 
04-Nov-2011 
16:00:00 - 


05-NOV-2011 
05:00:00-05:15:00 
-06:15:00 


DSN: Madrid, Spain Direction 

05-NOV-2011 


00:00:00-01:15:00 


Figure 1. Illustration of the orbit trajectory and the instances of measurement observations from the DSN ground 
stations and the TDRSS satellites during the first orbital period 
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For the second test case scenario in comparing the PF to the GMM, the same orbit is used. The 
measurements are available for the last 15 minutes of the orbital period toward apogee using the DSN in 
Canberra, Australia. In comparing the PF to the GMM, we are interested in looking at the state propagation 
of the PDF over time as it evolves through the orbit. The state PDF update performed at the final time 
will be illustrated in the comparison of the PDF capture representation between the PF and the GMM. 

VII. Simulation Results 


VII. A. EKF vs PF 

In Figure 2, we illustrate the estimation errors and the 2 -a error bounds for the EKF compared to the PF 
over 1 orbital period. The overall mean estimation errors (both prediction/propagation and measurement 
updates) are given in Table 1. The values in the table clearly show that the PF’s mean errors are closer to 
zero compared to the EKF. This can be attributed to the diversity of the particles that take into account 
the modeled uncertainties and the nonlinearities of the dynamic model. While for the EKF, the mean and 
the covariance only carry information up to the second moments alone and hence a lot of the uncertainties 
are pronounced during the prediction step but once the measurements are available, the EKF manages to 
converge the errors closer to zero, as is illustrated well in Figure 2(e). The values of the highest measurement 
update errors over one orbital period are shown in Table 2, for each of the filters during the measurement 
update instances. The errors in the tables and the plots are given in the VNB frame, which represents 
the transformed states from Earth-Centered-Inertial (ECI) to the body fixed VNB frame in the directional 
Velocity, Binormal and Normal components. The Velocity component is in the direction of the velocity 
vector, the Bi-normal component is perpendicular to the Velocity component and the Normal component is 
in the normal direction of the trajectory plane. 








PERIOD °' 6 


(a) PF: V-direction 


(b) PF: N-direction 


(c) PF: B-direction 



(d) EKF: V-direction 


(e) EKF: N-direction 


(f) EKF: B-direction 


Figure 2. Position errors (blue solid lines) and 2 -a error bounds (black dashed lines) for 1 period in VNB frame for 
PF (1000 particles) and EKF. 


Table 1. Mean position absolute errors in 
meters in the VNB direction over 1 period 


Direction 

Velocity 

Normal 

Bi-Normal 


EKF (m) 
3.075 
6.743 
2.834 


PF (m) 
0.0880 
0.631 
0.0758 
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Table 2. Highest measurement update 
errors in meters in the VNB direction 
during all four observations 


Direction EKF (m) PF (m) 
Velocity 6.532 6.768 

Normal -37.34 30.62 

Bi-Normal -19.31 —7.761 


We also investigate to see how the performances for the EKF and the PF compare to one another over an 
additional 4 periods of state propagation. In Figure 3, we illustrate the estimation errors and the 2-cr error 
bounds for the EKF compared to the PF over 1 orbital period and an additional 4 periods of propagation to 
illustrate the predictive accuracy of the PF. In this simulation, we see that for both the PF and the EKF, the 
position estimation errors are spiked at periapse in the Velocity direction. This is attributed to the fact that 
the dynamics of the orbital model are changing rapidly at this location and this is when the magnitude of 
the velocity is at its highest value. Moreover, we can also see the 2 -a error bounds for the Velocity direction 
are at their highest values at periapse and apoapse. A vice-versa behavior is noted in the Binormal direction, 
whereby the uncertainty bounds are at their lowest at apoapse and periapse, due to the law of conservation 
of angular momentum. Based on the fact that we are using 1000 particles for the PF, it is implied that the 
more particles we use the more accurate the PDF. So, we can assume that for an accurate prediction of the 
state PDF over longer periods of time, we can expect reduced growing errors especially around peripase and 
hence have a better way of predicting the state compared to using the EKF. The highest propagation error 
peak values are given in Table 3. It is also important to note that describing the errors bounded by 2-cr error 
bounds is not an accurate representation since the particles do not possess a Gaussian distribution, as will 
be illustrated in Sections VII. B and VII. C. 








Figure 3. Position errors (blue solid lines) and 2-cr error bounds (black dashed lines) for 5 periods in VNB frame for 
PF (1000 particles) and EKF. Note that apoapse occurs at day boundaries and periapses occur at half-day points. 


VII. B. GMM vs PF 

In Figure 4, we illustrate the state PDF estimate with measurement updates at the end of the first orbit. 
The GMM is illustrated using 1000 Monte Carlo (MC) runs to sample the state PDF and these samples 
are propagated through the dynamics to depict the Gaussian mixtures PDF capture. The Final PDF for 
the GMM is shown to have split from the initial single Gaussian PDF into three Gaussian mixtures with 
their respective weights given for the mean points GM1, GM 2 and GM 3 in Figure 4(a) and 4(d). The 
probability contours attempt to capture the probability values of each of the individual samples/data points. 
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Table 3. Highest propagation error peak 
values during the last four orbital peri- 
ods 


Direction EKF (m) PF (m) 
Velocity 269.9 186 

Normal 19.46 1.5 

Bi-Normal 8.998 4.988 


It is clear that the three Gaussian mixtures provides an approximation of the entire PDF using the Gaussian 
mixtures. In Figure 4(b) and 4(e), we illustrate the Gaussian mixtures using five Gaussian components with 
their respective weights given for the mean points GM1, GM 2, GM 3, GM 4 and GM 5. The splitting GMM 
has a library of the number of Gaussian components as well as the component weights, means and a fixed 
variance for each component provided in Reference. 24 One assumption is that the Gaussian components’ 
equal spacing from each other (hence the pre-calculated means of the components), suffices in capturing the 
non-Gaussian distribution. 

Moreover, a separate filter (the UKF in this case) is required to perform any measurement update, adding 
the extra computational cost for performing the UKF for each of the three Gaussian mixtures. On the other 
hand, the PF provides the state PDF as a collection of the individual 1000 particles with their respective 
updated weights using the information available from the measurements. Hence, each particle contributes 
to the capture of the PDF and there exists little to no overlap compared to the overlap that might appear 
using the GMM once we perform a measurement update. Since, the measurement noise is usually small, the 
updated PDF tends depict a more kurtosed distribution as illustrated by the PF, hence for measurement 
updates, the GMM and the UKF will provide Gaussian mixtures with extra overlap to account for the 
distribution present in the heavy tails. 

The colored particles are weighted Wi as follows: violet is for Wi < 0.002, cyan is for Wi < 0.004, blue 
is for Wi < 0.006 and black is for Wi > 0.006. In Figure 4(a), 4(b) and 4(c), the plots illustrate the PDF 
representation in the X-Y plane and in Figures 4(d), 4(e) and 4(f), the plots illustrate the PDF representation 
in the X-Z plane. 

VII. C. Non-Gaussian Measures 

To demonstrate the non-Gaussian nature of the PDF, we illustrate the higher order moments of the PDF 
propagation by showing the skewness and kurtosis of the state PDF over one orbital period in Figure 5. 
Skewness is defined as the third moment of the PDF, representing the shift of the PDF peak to the left or 
to the right of the center point. In Figure 5(a), we see how each position component’s skewness changes 
in the VNB direction. The Bi-Normal component shows more pronounced changes in skewness especially 
as it passes through periapse at the 48th data point. The Velocity component’s skewness shifts along the 
orbit trajectory from positive to negative as it goes from apoapse to periapse respectively. While, for the 
Normal component, we do not experience a profound pronounced skewness and this is due to the fact the 
orbital motion takes place in a plane hence there is no state evolution in the Normal direction apart from 
process noise. Kurtosis is defined as the peaked-ness of the PDF; a negative value represents a smoother 
peak and a positive value represents a sharper peak. This is computed as the fourth moment. It is clear that 
the Normal component is less peaked compared to a Gaussian distribution while the Binormal and Velocity 
components are highly peaked. In the Velocity component, we see that the distribution is highly kurtosed 
when the true anomaly is near 90° and near 270°. Also, in the Binormal direction, we see the distribution is 
highly kurtosed near periapse and apoapse relative to the rest of the locations on the orbit trajectory. For a 
Gaussian distribution, the skewness and kurtosis are zero and three respectively. In Figure 5(b), the excess 
Kurtosis (Kurtosis minus 3) is plotted. This represents the change from a Normal distribution, which would 
have a Kurtosis of 3. 


VIII. Conclusion 

Estimation of the PDF using the PF provides us with a more accurate representation of the region of 
uncertainty of our state estimates. For collision avoidance, this PDF knowledge will give us a more accurate 
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(a) 3-GMM: XY-plane 


(b) 5-GMM: XY-plane 


(c) PF: XY-plane 





(d) 3-GMM: XZ-plane (e) 5-GMM: XZ-plane (f) PF: XZ-plane 

Figure 4. PDF Position densities for GMM using 1000 MC runs and showing the 3 and 5 Gaussian Mixture’s mean 
location and contours compared to the PF’s particles and their weights to sum up the entire PDF. 
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(a) Skewness. 

Figure 5. The Skewness and Excess-Kurtosis of the 



(b) Excess-Kurtosis. 

over 1 orbital period using the PF for 1000 particles. 


representation of the likelihood of low-probablility but important events, that can be used with a decision 
making scheme to predict potential collisions and hence, implement the correct maneuvers required. The 
comparison of the PF to the EKF demonstrated that the PF has more potential and versatility for providing 
lower estimation errors and better propagation values with lower error bounds over time. In comparing the 
PF to the GMM, it was shown that the PF provides a more accurate PDF capture with less redundancy, 
compared to the GMM that uses a weighted sums of Gaussians that possess a potential of uncertainty 
regional overlap and redundancy. 
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